#!/usr/bin/perl -w
package Statistics;
use warnings;
use strict;

## user modules
use SimpleError;

sub is_numerical {
    $_[0] =~ /^
        [+-]?
        (?:
            \d+(?:\.\d*)?
            |
            \d*(?:\.\d+)
        )
        (?:
            [Ee][+-]?\d+
        )?
    $/x;
}

sub median {
    my ($ra, $sorted) = @_;
    error("not an array reference") unless ref $ra eq 'ARRAY';
    return undef if (@$ra == 0); # empty array
    unless ($sorted) {
        my @array = sort {; $a <=> $b; } @$ra;
        $ra = \@array;
    }
    my $med_ind = int(@$ra / 2);
    if (@$ra % 2 == 0) {
        # even number of elements in the array
        return ($ra->[$med_ind] + $ra->[$med_ind - 1]) / 2;
        #return [ $ra->[$med_ind], $ra->[$med_ind - 1] ];
    } else {
        # odd number of elements
        return $ra->[$med_ind];
    }
}

sub mean {
    my ($ra) = @_;
    error("not an array reference") unless ref $ra eq 'ARRAY';
    return undef if (@$ra == 0); # empty array
    my $sum = 0;
    for (@$ra) {
        $sum += $_;
    }
    return $sum / @$ra; # floating point division by default
}

sub rms {
    my ($ra) = @_;
    error("not an array reference") unless ref $ra eq 'ARRAY';
    return undef if (@$ra == 0); # empty array
    my $sum_square = 0;
    for (@$ra) {
        $sum_square += $_ * $_;
    }
    return sqrt($sum_square / scalar @$ra);
}

sub stdev_common {
    # a stdev is got as follows:
    # 1. calculate the mean of the population of numbers, which
    #    is actually regarded as the "expectation" of the numbers
    # 2. for each element of the number population, calculate
    #    the square of the difference of the number and the mean
    #    from step (1)
    # 3. calculate the sum of the numbers from step (2)
    # 4. calculate the square root of the sum from step (3),
    #    that's called "Standard Deviation"
    my ($ra, $mean) = @_;
    $mean = mean($ra) unless defined $mean; # step (1)
    my $var = 0;
    for (@$ra) {
       $var += ($_ - $mean) ** 2; # step (2) along with the summing
                                  # for late averaging
    }
    #return sqrt($var / @$ra);     # step (3), the mean along with
                                  # the square root calculation, or
                                  # step (4),
    return $var;
}

### standard deviation for an entire population
sub stdev {
    my ($ra, $mean) = @_;
    return undef unless (@$ra > 1);
    return sqrt( stdev_common($ra, $mean) / @$ra );
}

### standard deviation for a sample of a population
sub stdev_s {
    my ($ra, $mean) = @_;
    return undef unless (@$ra > 1);
    return sqrt( stdev_common($ra, $mean) / (@$ra - 1) );
}

1;
